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AN ATMOSPHERIC DENSITY MODEL FOR APPLICATION 


•IN ANALYTICAL SATELLITE THEORIES 

by 

Alan C. Mueller 


1.0 INTRODUCTION. 

If a fully analytical satellite theory which includes 
the drag perturbation is to be successful, it must possess 
three important qualities. First, the theory should be based 
on a canonical formulism whereby one can use the powerful 
tools provided by hamiltonian mechanics. Secondly, the model 
used to describe the forces acting on the satellite must not 
be so simplified that the theory becomes only a mathematical 
exercise. Lastly, the resulting theory must be concise so 
that the accuracy gained outweighs the extra computer costs 
required to reach that accuracy. 

Scheifele (reference 1) has developed an analytical sat- 
ellite theory based on the regular, canonical Poincare-Similar 
(PS(j)) element?. This is a very powerful set of elements which 
are in an extended phase space and have an independent variable 
which is similar to the true anomaly instead of time (refer- 
ences 2, 3 and 4). A very accurate and concise satellite theory 
has been developed to include the first order, short period 
and secular perturbations of an oblate central body. The drag 
theory has been built on top of the theory. 

The assumption in Scheifele 's theory is that the drag force 
is tangential to the orbit and proportional to the square of the 
velocity magnitude of the spacecraft. The constant of propor- 
tionality, which is a product of 'the density of the atmosphere, 
the ballistic number, and the drag coefficient, was not specified. 



Since the lifting force relative to the drag force and the 
inertial velocity of the atmosphere relative to the satellite 
velocity are small, the model used by Scheifele is adequate 
for giving the d-ireotion of the retarding force due to the 
atmosphere. Thus an important contribution to the analytical 
solution was made. The report (reference 1) is a concentrated 
effort to canonically transform the forces into the PS(j) space 
and also place them in a form suitable for solution. Therefore, 
the direction of the PS(j> canonical forces has been determined 
but their magnitude was not completely specified. Also, the 
tools of hamiltonian mechanics were employed to appropriately 
transform the forces and reduce the size of the equations. 

Numerical studies were conducted to confirm the accuracy 
of the resulting satellite theory. In the tests, both analyt- 
ical and numerical orbit predictors assumed that the density, 
ballistic number (weight over projected area) and drag coef- 
ficient were constants. Results showed that the analytical and 
numerical solutions matched extremely well, verifying that the 
transformation and quadrature solution were computed properly. 
Due to the unique character of the PS<{> system, the equations 
which . describe the motion are relatively simple and thus sat- 
isfy the first and third above mentioned qualities. 

However, for most satellites there can be extremely large 
changes in the density o'f the atmosphere along the orbit. Even 
for small eccentricities (e = 0.02) the density can vary by 
a factor of 100 

A study has been made (see section 2.0) to determine the 
errors that result from assuming an average constant density 
as compared to a density model such as that developed by Jacchia 
(reference 5). The comparisons point to the fact that the 
constant density model is adequate only for orbits of very small 
eccentricities. But in all the cases the analytical orbit pre- 
diction using a constant density model was much closer to the 
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numerical prediction which used a constant density than the nu- 
merical prediction obtained by employing the Jacchia model. 

This implies that it is the density model, not the analytical 
solution method, which restricts the accuracy of the analytical 
theory . 

Therefore the intent of this report is to develop an 
adequate density model and discuss the implications the model 
will have on the analytical drag theory. As in Scheifele’s 
theory the ballistic number and coefficient of drag will be 
assumed constant. 
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PRECEDING PAGE BLANK NOT PIT.MBD 

2.0 ' INVESTIGATION OF CONSTANT DENSITY MODEL 

If one assumes that the density of the upper atmosphere 
is a constant, then the drag theory and its cqr-respc^nding com- 
puter program described in reference 1 are essentially complete 
and could be made available to the user. The question is wheth- 
er or not this assumption results in solutions with acceptable 
errors.. To answer this question a series of numerical experi- 
ments have been' carried out. 

Since the density is so strongly dependent on height, 
several orbits over a wide range of eccentricities and semi- 
major axes were chosen by which to test the assumption. Choos- 
ing orbits over a wide range of the other orbital elements is 
not a neccessity fpr testing the assumption. 

Three numeTioally integrated solutions for the position 
of a satellite after a given time were determined for each of 
the orbit test cases.. Ail solutions include the perturbation 
due to the oblate mass of the earth, but the solutions differ 
in their drag model. The reference solution uses an extremely 
accurate but complex drag model by determining the density 
above the oblate earth with the model developed by Jacchia. 

A second numerical solution was found by assuming that the 

density is a constant . The constant density chosen was deter- 

% 

mined by using the Jacchia model to compute the density at the 
coordinates of the semi-latus rectum point of the initial orbit. 
Lastly, a solution was obtained by completely neglecting the 
drag force. The "NO DRAG" and "CONSTANT" density solutions 
were then compared to the reference Jacchia solution and the 
results displayed in table I for each of the test cases. Since 
drag so strongly perturbs the in-track position, the differ 
ences in the solutions are given by out-of-track and in-track 
position errors. Also note that the out -of -track error is 
shown in meters while in-track error is given in kilometers. 

By comparing the position differences resulting from the 
constant density model to the differences obtained by neglecting 
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drag, one has a relative measure of the constant density model 
assumption. For instance, in all the test cases with small 
eccentricities the CONSTANT solution always results in smaller 
errors than the NO DRAG solution. This is because the height 
of the satellite does not vary a great deal in the orbit and 
thus the vehicle does not observe a large change in the density. 
However, the cases in v/hich the eccentricity is somewhat larger, 
the satellite sees very large variations in the density and the 
CONSTANT solution is no better than neglecting drag completely. 

The conclusion is that the constant density model results 
in a reasonable solution only for very small eccentricities. 

Even under these tight restrictions, the solutions are only 
labeled "reasonable", not accurate. Other factors such as the 
diurnal variation of the density cause the constant density 
assumption to be crude even for circular orbits. For these 
reasons, this report will concern itself with the goal of de- 
veloping an accurate density model which may be inserted in 
the analytical theory. 
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TABLE I.- EVALUATION OF CONSTANT DENSITY MODEL 


a 

(km) 


h 

P 

(km) 


Position Differences 
Out of Track(m) / In Track(km) 


NO DRAG 


CONST 




411 ./ 16.32 

11 ./ 0.8 

417 ./ 16.48 

^ 36 ./ 1.6 

757 ./ 27.28 

368 ./ 12.24 

3815 ./ 129.28 

1664 . / 46. 08 

29 ./ 2.81 

4 ./ 0.14 

29 ./ 2.84 

3 ./ 0.26 

43 ./ 3.88 

10 ./ 1.35 

92 ./ 7.68 

54 ./ 4.72 





























-17- 


IjlO DENSITY MODELS fjttijlJEDiNG PAGE BI4ANK NOT FILMED 

In developing a density model for the analytical theory, 
one is severely restricted by the fact that the model must be 
in the form of a fourier series in the true longitude. As is 
the case in most analytical theories, the perturbation must be 
written in a fourier series to facilitate the solution of the 
differential equations of motion. Usually, the solution is 
obtained by quadrature. 

Several density models have been developed to predict 
accurately the density at any point in space and time. Examples 
are the Jacchia model (reference 5) and the USSR model (refer- 
ence 6). But both models are extremely complicated and too 
unwieldly for analytical satellite theories. In the analytical 
theory of Brouwer and Hori (reference 7) , the density was 
assumed to be an exponential function of the position radius 
of the satellite. The exponential must then be expanded in a 
Poisson series so the quadrature can be performed. This model 
has several difficulties. It first has a problem with conver- 
gence, which Brouwer points out. Secondly, it is simply a poor 
model for describing the dynamic atmosphere. The density is 
extremely effected by such factors as the level of solar activ- 
ity and whether it is summer or winter, day or night. Thus the 
model in Brouwer, Hori theory is simply inadequate. 

Recently an extremely simple density model (referred to 
here as B-M) has been developed to match the Jacchia model to 
a high degree of accuracy (reference 8). The variations in 
the density due to changes in the height and changes in the 
relationship of the vehicle and sun position (diurnal effect) 
are included explicitly. Long period variations such as the 
changes in the average solar activity and semi-annual vari- 
ations are included implicitly in the coefficients of the model. 
The value of the fcoef f icients are determined by a procedure 
called "calibration”. The simple formulation allows the model 
to be inverted, i.e. given the density at different points in 
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space (as determined from Jacchia) one can compute the coeffi- 
cients ^ of' tlie^B-M model'^.''* * Since the coefficients are implicit 
functions of long period effects, they can be considered 
constants over a limited period of about a month. 


Even this extremely simple model cannot be applied in 
the analytical drag /theory because it cannot be written in the 
form of a fourier series. However the technique of the B-M 
model does give important insight and direction to follow. 


In all the models discussed, the representations of the 
kensity are considered to be global. In other words, given 
any position in the atmosphere one can determine its density. 
The approach, to be taken here, is to develop a model which 
expresses the density along a particular orbit. The coeffi- 
cients in the model will be calibrated with the Jacchia model 
in a manner similar to the B-M model. But in this case, the 
coefficients in the model are not only implicit functions of 
the long period variations in the atmosphere but also the 
orbital elements which describe the orbit. The result is a 
density model which can be written in a fourier series and 
easily implemented into the drag theory. Since the orbital 
elements are perturbed by drag, the coefficients in 

the model must be updated periodically or corrected in some 
manner to reflect. the changes. Since the perturbations are 
small, the updating would be infrequent. 


3.1 Development of the Model 

In the B-M model, the variations in the density due to 
the height and diurnal effect are modeled as 

(tJH) -I- (1) 

diurnal terra 


p = exp 
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where 

p = density 

T'(h) = function of altitude (nighttime vertical profile) 
S(h) = function of altitude (diurnal magnitude) 
g(as , 6^) = function of right ascension and decli- 
. nation of the sun and vehicle. 

And similarly the USSR model is expressed as 

p = exp (T’(h)) * (1 + S'(h) * g"" (a ,6^,a .6^ ))• (2) 

S S V M 


where 

T'(h) = function of altitude (nighttime vertical profile) 

S'(h) = function of altitude (diurnal magnitude). 

Both models point to the fact that the density may be expressed 
in such a form as 

p = T*(h) + S*(h) • g"" (a^,6^,a^,6^). (3) 

where T* and S* are functions of the height and g is a 
simple function of the angular coordinates of the sun and vehi- 
cle. If the oblate figure of the earth is neglected, then T* 
and -S* may be assumed to be functions of the vehicle position 
radius. 

i* 

in the PS(|) ' element system the radius is described as 
follows 


P 

r - 

1+e cos (j) 



(4) 


where p is approximately the semi-latus rectum and n is 
proportional to the eccentricity, and is given by 

+ 

The PS(i) action variables are > p^ , p^ and their 

canonical conjugate variables are 0^, 0.^ » 0^- 
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= $(p2 eos sin a^) 

The radius can be expanded about p 
considered small, i.e. 


3 ^ 


V 


if the eccentricity is 


n 

r = p ®i^^^ ^1 

i=o 

where e^ are functions of p 

Since T* and S* are functions of the radius, this 
suggests they too may be described in a similar manner 

n 

T*(h) =2 a' e 
i=o 

n 

S*(h) = Z) t>’ cj 
i=o 

where the coefficients a.^ and bl are implicit functions 
of the eccentricity, semi-latus rectum p , and the character 
of the atmosphere . 

Neglecting small terms in the angular function 
g(a ,6 , oi ,6 ) the USSR model gives 

s S V V 

m 

( 1 + <?0 8 \ 2 

-T-) 

where 

Gosip = ^ (z • sin 6 + Gos 6 • (x • oos y + y * sin y)) 

X s s 

y = + (p 

x,y ,2 = defines position of satellite 
(a^,6^) = right ascension and declination of sun 


( 6 ) 

( 7 ) 
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(j) = defines lag of the density bulge behind the 
sun C(})=^37®). 

The power of exponent m ranges near the value of 4. 

If m = 4 is chosen then g“(a) can be expressed as a fourier 
series in the PS<^ elements 

d» = g**(a) = dQ + d^ sin oos 

+ d^ sin 2a^ + d^ cos 2o^ (9) 


where 


, 2+D^+E^ 

% = 8 


D 


E 

*4 


d 


l' 


E 

2 


(D^-E^) 
4 8 



D = + oos 6 oos Y 

3 . s 

E = - S-B + oos 6 sin y 
3 s 


( 10 ) 


B = sin 6 /2(G+H)‘ 
^{j L s 


- cos Y P3 


sin 



The coefficients d^^^ can be considered constants over a few 
revolutions'. But due to the fact that the position of the sun 
changes and the orbital elements are perturbed, an infrequent 
update will be required to reflect these changes. 

Thus the total model for the density along the initial 
orbit follows from equations (3), (O), (7) and (9) 

n 

Pq = S (a. +' d b ) ? 5 ; 
i=o 


( 11 ) 
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where 


a. = a! + d^bl 

X X ox 

d = d' - d 



( 12 ) 


The coefficients a. and b. can be determined by a cali- 

XX 

bration with Jacchia. A linear system of 2(n+l) equations 
with 2(n+l) unknowns must be solved to determine the coef- 
ficients . 


3.2 Density Model Corrections 

The density model proposed in section 3.1 reflects the 
observed density variations due to the two-body changes in the 
height and due to the varying, angle between the sun and the 

f 

satellite. Santora (reference 9) also considers two additional 
phenomena that effect the height of the satellite which in turn 
causes variations in the observed density. One effect is th^ 
changes in the height due to oblate figure of the earth. For 
instance, a satellite in a circular orbit about the equator 
will see no changes in the height, but a satellite in a cir- 
cular, polar orbit will find that its height will change by 
about 20 km. These changes in height translate to a substan- 
tial variation in the observed density of about 50% 

Secondly, changes in the height due to the periodic 

oscillations in the radius can also result in large variations 
in the observed density. 

Both of these effects may be modeled as corrections to 
the model proposed in section 3.1: 



v^’here Ah is the change in the height due to the above men- 
tioned effects, is given by equation (11) and p is the 
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corrected density. The constant a is an implicit function 
of the character of the atmosphere and may be determined easily 
by calibration with the Jacchia model. For small changes in 
height, one can expand and truncate equation (13) which results 
in 


p = cPq , c = 1 + aAh (14) 

The oblate figure of the earth can be described by 

R = rJ (l+6(|)^)"^ (15) 

where R is the radius at an arbitrary point on the earth 
R^ is the mean radius at the equator 

6 is the bulge parameter 
z 

(— ) determine the latitude of the point . 

Since 6 is small (6- 0.67x10 ^) equation (15) can be ex- 
panded and truncated to give 

R = \ (16) 


-The latitude term may be expressed in the PS(|) 
G+H / 

* <=3 203P3 Bin 

-I- (o^-pp aos 


element as 


(17) 


Thus if one defines the mean radius as observed by a satellite 
in its orbit as R 

m 


R 

m 



(G+H) 

8G^ 




1-6 


(18) 
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Then the oscillation of the height about the mean radius 
is given by 


Ah^ = R 
R € 


6(G+H) 

8G^ 




203P3 


sin 2cr, + 




cos 2a. 


( 19 ) 


Derivation of the change in height due to the 3^ term 

is somewhat more lengthy but can be simplified by neglecting 
terms on the order of the eccentricity. Since the oscillations 
in the height due to have a magnitude of about 5 km and 

•becg,use the drag theory is restricted to e== 0.1 , neglecting 

order 0(e) terms results in errors of at most 500 meters. 

This results in a negligible error in the density computation. 

The change in height due to can be found by differ- 

encing , the radius computed from the mean elements r' from 
the actual radius r : 

Ar = r - r’ (20) 

where 


P' 

r' = 

1+Q'Z' 


P 

r = — = 

1+QZj 


Primed variables are based on the primed (mean) PScj) elements 
(see reference 4). Neglecting 0(e) terms, the difference 
■becomes 


Ar = p(l-QZp - p'(l-Q'Zp (21) 

If one defines relations between primed and unprimed values as 

p = p ' + Ap 
Q = Q’ + AQ 


( 22 ) 
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then Ap, AQ and can be found from relations derived from 

the generating function used in eliminating short periodic 

effects (reference 4). 

From the relation for the semi-latus rectum p in terms of 
PS(J) elements, 

(23) 

one can linearly approximate Ap as 

+ PjAPj + — — j- Ap^ (24) 

(2P4)2 

The deviations Ap^ and Ap^ can be found directly from 

the partials of . But Ap is of the order 0(e) since 

and p^ are eccentricity terms and Ap^ = 0 because of 

the use of total energy elements. A similar argument can be 
made for AQ . And thus Ar reduces to 

Ar = pQAZ^ + 0(e) (25) 

From the expression for 

Z^ = p^ oos sin (26) 

one can approximate AZ^^ as 

AZj^ = Ap^ oos - Ac 2 sin 

- ACj^ (P 2 sin oos ct^^) 




(27) 
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I^eglecting terms of 0(e) 

AZj^ = Ap^ Gos - Ao 2 sin 


( 28 ) 


By a rather lengthy but straight forward derivation one finds, 
through use of the generating function , that AZj^ can 

be expressed as 

w 

AZ, = e (G^- 3H^ + s sin 2a, + c oos 2a.) (29) 

1 3G2 V 1 1/ 


where e is the 
erence 4 and 

w 

s 

c 


perturbing parameter defined in ref- 


= _Q_ 

2pq 

= - (G+H) a^pg 


= (G+H) 


2 


(30) 


If one defines the average radius r^ as 

epQw 

r = r' + (G^-3H2) (31) 

^ 3G^ 

then the oscillation in the height due to J 2 is given by 
-epQw 

Ah = (G+H) (2a_p sin 2a, + (a^-pM 00 s 2a, ) , (32) 

^2 qq2 \ j j i 33 1/ 

It is interesting to note the similarities in the equation 
for the change in height due to the dynamics (eq.(32)) and the 
geometry (eq.(19)). Both the oscillations will have two peaks 
and two valleys, except their magnitudes are different and 
are 90® out of phase. Also, both oscillations vanish when the 



~2'7- 


inclination goes to zero. Thus the model given by equation 
(14) can be expressed as 


1 + 


K ^20^9^ sin 2 a ^ + (a^-Pp oos 


( 33 ) 


whQre 


(G+H) R 6 epQw 

K = a 

2G^ 4 3 


(34) 


The coefficients in Pq (equation (11)) are computed with the 

primed PS(}) elements by a calibration' to the Jacchia model, 

as described in Section 3.1. The height used to evaluate the 

Jacchia model is given by the difference between the averagfj 

radius r and the mean radius of the earth seen by the sat- 
a 

ellite 


h = 


r 

a 



epQw 

h = r' + (G^-SH^) 

3G^ 


R 


6 (G+H) 

1 p (a^tpp 


8G" 


L 


(35) 


3.3 Model Verification 

To verify that the proposed density model has an adequate 
accuracy, a set of experiments have been conducted to compare 
the new model with the Jacchia model. 

In all of the tests, the density as observed along dif- 
ferent points of a pei*turbed orbit about the oblate earth 

is computed with the Jacchia model and then compared with the 
new model. A value ’ of 4 was chosen for n (eq.(ll)) so that 
pQ has 10 coefficients to be determined by calibration. A 

matrix inversion algorithm (reference 10) was used to solve for 
the value of the coefficients. 
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Several different orbits and sun conditions were chosen 
by which to compare the new models. These cases have alT been 
labeled in Table II and will be referred to in the discussion 
of the comparisons . 

The first comparison of the two models is intended to 
demonstrate that extremely large variations in the density 
due to the solar activity level, ^ , can be accounted for 

by the new model's technique of calibration. 

In figure 1(a), the density is plotted versus the true 
longitude of the satellite. Two graphs are shown, one as 
determined from Jacchia, and the other as found from the new 
model. The initial conditions are given by Case 1 in Table II. 

I 

Jhe resulting differences between the two models are never 
greater than three percent. In figure 1(b), another two 
graphs are shown except with initial conditions given by 
Case 2. Again the models are in good agreement. The only 
difference between case 1 and 2 is their solar flux values; 
but notice the large quantitative differences between the 
density plots in figures 1(a) and 1(b). The new model is able 
to account for these very important differences. 

The second comparison demonstrates the ability of the new 
model to account for density variations due to the changing 
position of the sun. Again in figure 2, the density is plotted 
versus the true longitude. There are four plots on figure 2 
representing the two different models and two different initial 
conditions given by case 3 and 4, The only difference between 
these cases is the position of the sun. In case 3 the siin is 
in the first day of summer while in case 4 the sun is in the 
first day of spring (vernal equinox). Here too one finds that 
the new model is able to reflect the differences. • 

The four remaining comparisons are intended to show that 
the proposed model simulates the Jacchia model for a variety 
of orbits and to demonstrate the errors resulting from neglect- 
ing the effects of the diurnal bulge, the oblate figure of the 
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earth, or the periodic oscillations in the radius. In 

figures 3(a) , 4(a) , 5(a) and 6(a), the density as determined from 
Jacchia is plotted versus the true longitude, corresponding 
to the test cases 4 through 7. These peculiar and varied 
graphs give a clear picture of the problem of developing the 
density (vertical axis) as a fourier expansion of the true 
longitude (horizontal axis). Be careful to note the changes 
in the scale for the density axis for the four different plots. 

In figures3(b) , 4(b) , 5(b) and 6(b), the percentage difference 
between the new model and Jacchia is plotted versus the true 
longitude- In each figure there are four plots labeled A, 

B, C and D. The D plot is the differences resulting from 
the new model if the diurnal effect is neglected. Similarly 
the C plot is from neglecting the oblate figure of the earth 
and the B plot is from neglecting the radius oscilla- 

tions. Finally the A plot is the resulting differences if 
the complete model is used. 

Case 4 is a circular, polar orbit and thus the density 
variations seen in figure 3(a) are due entirely to the diurnal 
effect and the changes in the height because of the oblate 
figure and mass of the earth. The rather strange appearance 
of the plot is a result of the cancellation and addition of 
the different effects. From figure 3(b) one finds large errors 
result if diurnal or oblate figure effects are neglected. 

Smaller errors result from neglecting the oscillations. 

As expected, the oblate mass and oblate figure effects have 
two peaks and valleys and are 90® out of phase. Also note 
that the errors in plot C are always positive. The reason 
for this behavior is because the earth's radius was assumed 
to be a constant equatorial radius which is larger than at 
the poles. The result is that a greater density is predicted 
for every point except when the satellite crosses the equator. 

The residual errors from the total model are always less 
than three percent. 
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Case 5 is also a circular, polar orbit but has an. altitude 
of about 300 km. greater than in case 4. The density plot of 
this case is shown in figure 4(a) and the differences are plot- 
ted in figure 4(b). 

The neglecting of the diurnal effect results in the larg- 
est errors. The differences between the magnitude of the errors 
in plot C in figures 3(b) and 4(b) reflect the fact that the 
diurnal effect becomes more prominent for higher altitudes. 

The other plots show the same pattern as in case 4 and again 
the total model exhibits small errors. 

Case 6 is a polar orbit with what appears to be a rela- 
tively small eccentricity. Actually this small eccentricity 
translates, to a very large density variation. In figure 5(a) 
the major variation in the density is due to the variation -in 
height of the elliptical orbit. But by examination of the 
error plots in figure 5(b), one finds the other effects are 
also very important - especially the diurnal effect , 

Case 7 is also a polar orbit with an eccentricity which 
results in density variations of 4.22 x 10“^^ to 

6.03 X 10”^^ ^ , a factor of over 100. Here the plot in 

m^ 

figure 6(a) begins to look more like an impulse where the 
satellite experiences large densities only near its perigee 
point. The scale of the vertical axis in figure 6(b) has been 
changed to show the very large errors resulting from neglect- 
ing the diurnal effect . Plot B has not been shown because 
the scale does not bring out the differences between plots 
A and B 

Errors in plot A become as large as ten per cent but 
this is only when the satellite sees its smallest densities 
at the apogee . 

As the eccentricity of the orbit increases, the plot of 
the density versus true longitude begins to look more and more 
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like an impulse. The result. is that for eccentricities of 
e > 0,03 the model begins to break down. The fourier series 
cannot converge when modeling the impulse. 

However, this is 'not such a severe restriction. Most 
near earth satellites have eccentricities well within the 
^ = 0.03 . For satellites with greater eccentricities one 

may avert this problem in a manner similar to that described by 
Watson (reference 12) . The density can be modeled only in the 
region near the perigee point where it is largest . The ana- 
lytical integration then proceeds in steps. 

From these results, it is concluded that the proposed 
model is quite suitable for applications in the analytical 
theory. It matches the Jacchia model extremely well and is 
easily written in the form of a fourier series. 


3.4 Second Order Corrections 
— — , , ■ 

Because the density is so stro,ngly dependent on the height, 
it is important that changes in the height are accounted for. 

The proposed model accounts for the two-body and changes 

in height and also the changes due to the oblate figure of the 
earth. But the drag force itself causes a secular perturbation 
in the height. For very high satellites, the ' perturbation is 
small and can be neglected in the density considerations. But 
for satellites which pass through the dense atmosphere, the 
effect is much like a snowball growing as it falls down the 
hill. The drag force causes the satellite to dip deeper in 
the atmosphere where the more- dense air causes a stronger drag 
force and so on. Once again this effect can be modeled as a 
correction as in equation (13) 

p = (1 + aAh) • 


(36) 



- 39 - 


If one assumes that the drag perturbs the radius of the satel- 
lite orbit in a secular manner and neglects the terras on ,the order 
of the. eccentricity, then the change in height can be related 
to the perturbation in the semi-latus rectum, i.e. 


Ah = Ap 


( 37 ) 


Since the p^ is given by 

1 


P = 


~ iCcf^+Pp + 


u 


(2P^) 


then the deviations due to drag can be written as 



/ajAo^ + PjAPj + 


y 

(SP^)' 



( 38 ) 


But neglecting eccentricity terms and assuming a secular drag 
perturbation in the energy one finds 


Ap 




(2p^)’^2 At) 


T 


( 39 ) 


where 



is the rate of change in the energy due to drag. 


Initially 



may be determined by neglecting the effect 


At 

of drag on the atmosphere density model. Once this guess is 
made, one can re-evaluate the changes in the elements based 
on the complete model. 
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3.5 Density Model Summary 

Prom the theoretical developments in this section, the 
complete proposed density model may be written as 


p = 1 + ( a — 1 T + aAh p 


Ah = 

K j^2a^ 

sin 

2ai ■ 

Ap 

-1 

/yp ' 


At 

2P4 V 

2P, 

At 


G+H 1 

('K 6 
e 

epQw 

Is — 

2G2 \ 

L 4 

3 


-2 2 
f ^ — P „ 

3 




d = d^ 6-in + d 2 oos + d^ sin 2o^ + d^ eos 


E 

D • - 
4 


(D^-E^) 


D = CT„B + QOS 6 QOS Y 
3 s ' 
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E = -p„B + cos 6 sin Y 
3 s 

1 

B = 

2G 

Y = ttg + (|) 

a , 6 : right ascension and declination of sun 

s s 

: defines lag of diurnal bulge behind sun 
6 : measure of geometric oblateness of earth 
e ; perturbing parameter 

j : model parameters found by calibration 




sin 5 . \/2(G+H) - QOB 5 (a^ eos y + Po Y 

S S O ' o 
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4.0 DENSITY MODEL IMPLICATIONS ON DRAG THEORY 

The proposed density model requires a number of fourier 
expansions in addition to those developed for the constant 
density case in reference 1. The addition of the new density 
model in the analytical theory does not, however, require any 
new Taylor series expansions which were processed by the com- 
puter logic described in that reference. Thus the discussion 
here will be concerned with the new fourier expansions that 
are needed. 

In the discussion of the power series expansion, one must 
make a distinction between the expansions of the properties of 
the orbit (such as the radius, time and velocity) and the ex- 
pansion in the density model. Even though both sets are power 
series expansions in the eccentricity, the density model does 
not converge as rapidly as the orbit properties. Thus for 
small eccentricities one may truncate at a very low order for 
the orbit properties, but still be required to carry out the 
density expansion to a higher order. In addition, it is not 
clear, at this time, to what order one must carry out the prod- 
uct of these two different expansions. 

Postponing arguments on when and how to truncate, 
one may still list all the fourier series that are needed for 
a specific order m in the orbit properties expansion and n 
in the density expansion. They include: 

, cdc^ 

c?^ eos , cdC^ cos 

c?^ sin , cd^^ sin 

2 ’ 


(40) 
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s 

®"^” *^1 ’ ^1 

for 1=0 , l...n+m ^ where c is the corrective term 

given by equation (14) and d is the diurnal term (eq.(12)). 

^2 is given by the relation 

?2 “ ^ ( p 2 *^1 '*’ ^2 

If n and m are both set to 4 , then there are a 

total of 108 expansions needed. This seems at first to be 
a rather imposing number of expansions. But drag is locally 
a very small force and thus the periodic effects of drag may 
be neglected. Therefore, only the average of these expansions 
needs be known. The average of a fourier expansion is simply 

the zeroth order term and thus only that term must be computed. 
One exception to this rule is the evaluation of the mixed sec- 
ular forms in the time_element . In this case the full expan- 
i i 

sions for cC^ and cdi;^ must be evaluated. Therefore, only 
18 full expansions need be determined. 

If the following notations are made 


i+2 



( 41 ) 
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and if the mean of a function f(a^) with respect to a 
is defined by 






( 42 ) 


then the average values of the expansions of equation (40) 
are given by 





The average of the expansions involving the product of 

c and d , can be found by replacing the uncapped 

i i i 

Xj with the capped values ' 
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4.1 Determination of Fourier Coefficients 

The coefficients of the fourier series given by equa- 
tion (41) must be evaluated in some manner to complete the 
solution. There are several options one might take to accom- 
plish this t 9 .sk. This section will be devoted to exploring 
these different options. 

Perhaps the simplest option is to truncate the expansions 
on a small order of the eccentricity. This reduces the number 
of fourier series to be evaluated but also restricts somewhat 
the generality and accuracy of the solution. However, one 
should remember that most drag perturbed satellites do in fact 

V 

have a very small eccentricity. Also a solution with .a small 
error is better than no solution at all. If one reduces the 
equations by truncating the solution, then one may determine 
the fourier coefficients by explicit equations which are either 
derived by manual or computer manipulation. These explicit 
equations could be programmed without requiring an extremely 
large storage allocation. 

A second option is available whereby one can evaluate all 
the fourier expansions without the loss of generality or accu- 
racy. To eliminate the computer storage required, of the ex- 
plicit expressions, one can compute the coefficients of the 
fourier series in a recursive manner. The root of this proce- 
dure is an algorithm which, given the numerical value of the 
coefficients of two fourier series, can evaluate the coeffi- 
cients for the product of the two series. The algorithm is 
developed from several trignometric identities and is outlined 
in Appendix I . With such an algorithm, all the coefficients 
in equation (41) can then be found in a recursive manner. For 
instance, the coefficients for can be found from the co- 
efficients of and the coefficients for can be eval- 
uated from the coefficients of and , and so on. 

This reduces much of the instruction computer storage required 
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by explicit expressions. Besides the advantage of being com- 
putationally simpler, this method also allows for easy updating 
of the drag model to include terms which effect only the magn-t- 
tude of the perturbing force. Thus the density model can be 
changed easily to an expansion of arbitrary order. ■ More im- 
portantly, if modifications of the density model, ballistic 
number, or coefficient of drag can be expressed as a fourier 
series (in a^) , then the recursive algorithms allow for a 

much simpler implementation of these changes. The laborious 
manual derivation of these terms can be eliminated and thus a 
typical problem of analytical theories is partially avoided. 

The disadvantage to this method is the fact that no explicit 
equations are developed (which is ironic because this is what 
the method was designed to avoid) . The algorithm is numerical 
in nature and like all numerical solutions one can find the 
"correct" solution but loses the insight gained from analytic, 
explicit solutions. In addition, the mathematical solution 
and the computer algorithm are so interwoven that they become 
inseparable. This leads us to the third option. 

A recursive solution for computation of the fourier series 
may be, possible without- the loss of insight encountered with 
option two, above. In Mueller (reference 11) the perturbing 
geopotential is written in the form of a fourier series by the 
use of recursive relations. By a similar approach, the per- 
turbing drag canonical forces may be developed in a recursive 
manner using literal expressions throughout. Thus the solu- 
tion would not have the "loss of insight" typical of numerical 
methods and would not be tied to a computer algorithm. As an 
illustration, a recursive expression for the powers of 

has been developed in Appendix II. Since the powers of 

are the major components of the fourier series given by equa- 
tion (40), the expressions of Appendix II more clearly define 
and justify the approach. 
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The fourth option, and the one the author feels is most 
viable, is an approach which is a careful blend of all the 
options discussed. A "careful blend" would be one that brought 
out the advantages of each of the options and minimizes the 
disadvantages . 


4.2 Qualitative Aspects of Drag Perturbation 

By examining the canonical PS4> forces in the light of 
the new density model, one may ascertain a qualitative descrip- 
tion of the effects of drag on the elements. 

The PS^ elements (the energy), p^ and 

(related to the inclination) are perturbed secularly by the 
static density profile on an order of 0(y) . The static den- 

sity profile refers specifically to that part of the density 
model which contains the a^ coefficients (eq. (11)). y is 

the ratio of the magnitude of the drag force to the two-body 
force magnitude- The eccentricity terms, and p^ , are 

secularly perturbed on the order of 0(yc) where e is eccen- 
tricity. The diurnal terms (the part containing b^ coeffi- 
cients) cause additional secular perturbations of 0(ye) in 
all the elements (except which is not effected at all by 

the drag force ). Also, the static density and diurnal terms 
(Pq in eq. (11)), give rise to quadratic and mixed secular 

terms of 0(y) in the time element . It is important to 

note that the neglection of the diurnal term results in errors 
of 0(ye) in all the elements except . The diurnal term 


In reference 1 it is shown tha.t the motion of (related 

to the true longitude) is not affected by drag perturbations. 
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dooK no ( averai^o out in the along track position error 

The correction for the static- profile for the earth's 
oblatonoss and perturbations of the satellite's altitude 

la'sulL in so<uilar perturbations OCcye^) where c is the 
magnitude of the correction (c.10.5) . The reason for the 

iu^noro of the ec.cen tr i c i ty is th(' fac.t that the correction is 
a sinusoidal function of 2o ^ . But what is interesting is 

I he coupling of the' correction and the diurnal terms whic^h 
results in terms of ord<>r 0(cy) . Thus the correction does 

» I 

not vanish for vanishing eccentricities. 

The correction of the density model due to drag results 
in quadratic and mixed secular perturbations of order O(c-y^) 
in all Hie elements except . The mixed se<iular term.s 

may be neglected because of the size of the perturbation. In 
the time element , the correction gives rise to c;ubic and 

mi.xed quadratic terms of order O(cy^) 
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fREOEDlNG 

5.0 NUMERICAL STUDIES 

Most, but not all, of the density model developed in 
section 3 has been implemented in a prototype program for the 
analytical drag theory. This new program is a revision of the 
program developed in reference 1 for the constant density mod- 
el. It resides on the UNIVAC 1110 and may be executed in 
the interactive mode without the necessity of an overlay. The 
second method described in section 4.1 and outlined in Appen- 
dix I was used to generate all the necessary fourier coefficients. 

Because time did not allow, the corrections in the density 
model due to the earth's oblateness and perturbations were 

not included in this program. But to test the remaining parts 
of the drag model, only equatorial orbits were chosen for test 
cases. For equatorial orbits, the neglected terms become very 
small and thus the errors that are observed in the tests should 
be realistic. As in the numerical studies of section 2, an 
extremely accurate numerical solution was generated which in- 
cludes the oblateness and a precision drag perturbation. 

The density used in computing the drag force was obtained from 
the Jacchia model. The coefficient of drag’ was set to C^=2.2 

and the ballistic number is an average value for the shuttle 

B = 100 ib/ft^ This numerical reference solution was then 

n 

compared to three different analytical solutions, all of which 
include the short period effects of - One solution ne- 

glected drag completely, (NO DRAG) , while a second solution 
included drag but considered the density as a constant (CONST) 
The third solution was computed with the new program with the 
proposed density model (TOTAL) . The results are tabulated 
in tables III and IV. The size and shape of the different 
orbit test cases are given by the semi-major axis a and the 
eccentricity e (h^ is the perigee altitude). The position 

errors between the analytical and numerical solutions are 
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displayed by the out of track error (given in meters) and 
the in-track error (given in kilometers). The solutions were 
compared about one half day or eight revolutions from the epoch 
of initialisation!. 

Results of satellite orbit prediction experiments are dis- 
played in table III. The first case is a very low, circular 
orbit where the drag force is extremely large and the diurnal 
effect is small. Therefore the major difference between the 
CONST and TOTAL solutions is the method in which the density 
model is corrected to account for the drag's lowering of the 
altitude. This correction does give a much better solution. 

The next four cases, with initial semirmajor axis of 6678 km 
shpw that the TOTAL solution gives a general improvement 
over the NO DRAG solution of better than one digit in the 
out of track position and almost two digits in the in-track 
position-. The improvement over the CONST model is substan- 
tial for the in-track position but is not so great in the out 
,of track range.. In .fact, for circular orbits the CONST so- 
lution appears to be better than the TOTAL model. This in- 
consistency is most probably due to two interacting effects. 
First the constant density model results in out of track errors 
which are on the order of the eccentricity. This explains why 
the CONST and TOTAL errors are of the same order for small 
eccentricities. The long period effects of are also of 

the same magnitude for this relatively high satellite. These 
long period effects have been neglected in the analytical solu- 
tions, and thus they corrupt the -estimates of the errors in 
the drag model. For satellites which have even higher alti- 
tudes, the o,ut 'Of track error is due almost entirely to the 
long period effects of . For -this reason, the test cases 

with a larger semi -major axis are not displayed. One can infer 
from these empirical results that Ihe drag model has been re- 
fined to such an accuracy that long period effects should 

be included to maintain a consistent accuracy. 



-53- 


Another numerical test was carried out to demonstrate 
that the diurnal effect is indeed modeled correctly. An orbit 
with a semi-major axis of a = 6878 km and eccentricity of 
e = 0.04 is used as the test case. This represents an orbit 
with a perigee height of 224 km and apogee height of 775 km. 

The relation of the orbit with respect to the sun is a strong 
factor in determining how much the satellite will be perturbed. 
Again three analytical solutions are compared to a reference 
numerical solution. However, for these tests none of the so- 
lutions include the perturbation, so that the errors are 

purely from the drag model. The results are shovm in table IV. 
The Jacchia model was used for the reference solution and 
assumes that the sun is at the vernal equinox. Two TOTAL 
solutions are computed; one with the sun at the vernal equinox 
(SPRING) and another which assumes the sun is at the first 
day of summer (SUMMER) . As one can see the TOTAL (SPRING) 
solution is certainly the most accurate. A large error is 
incurred if the sun is assumed to be in summer when actually 
it is spring; as seen from the results of TOTAL (SUMMER) 

The numbers in parentheses are the errors when the TOTAL (SUM- 
MER) solution is compared to a reference solution where the 
Jacchia also considers the sun to be in the summer. 
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TABLE III.- NUMERICAL SOLUTION vs ANALYTICAL SOLUTION 


r 

a 

(km) 

e 

h 

p 

(km) 

Out of tra< 

3 k (m) / In t 

/ ' ^ 

rack (km) 

NO DRAG 

CONST 

TOTAL 

6578 

O.Q 

200 

4891, /172. 8 

710./12.0 

15.2/0.3 

6678 

o:o 

300 

411./ 16.32 

11./0.83 

40.6/0.39 

0.001 

293 

417./ 16.48 

36. /I. 63 


0.01 

233 

757./ 27.28 

368./12.24 

50.0/0.008 

0.02 ! 

l^H 





TABLE IV.- DIURNAL EFFECT 


Position Error 

NO PRAG 

J 

TOTAL 

CONST 

SPRING 

SUMMER 

Out of Track 

3279. 

169. 

15. 

197. (11.) 

(m) 





Iq Track 

59.76 

25.12 

2.32 

13.6 (1.68) 

(km) 












































6.0. CONCLUSIONS 


In the past-, the computation of the drag perturbations 
by an analytical method has not. been feasible for.- two reasons. 
First, the various element sets chosen to base the past drag 
theories have never resulted in a tractable set of differential 
equations. This is a requirement in order to obtain • concise 
expressions for the solution. Secondly, the analytical drag 
theories have failed because of their inability to model the 
"real world" density as a fourier series. Both of these prob- 
lems are eliminated through the development of the drag dif- 
ferential equations in the PSij) elements, and the subsequent 
development of a new density model, which closely matches a 
realistic density model. 

The new density model has several distinct advantages. 

It has a rather subtle yet extremely important advantage in 
that the model is a power expansion in : i.e., the basic 

expansion used in Scheifele's development of the differential 
equations. The density model is much easier to implement in 
Scheifele's theory than, for instance, a model developed in 
power expansions of the radius . That type of model would 
require additional equations in order to be placed in the 
basic form of the expansion in powers of 

Another advantage of the new density model is its abilit 3 ? 
to simulate density model. In this report the Jacchia 
model was chosen for the simulations, but any other model 
could have been used. If more accurate density models are 
developed, they may also be chosen for simulation and thus 
the analytical drag theory can reflect the accuracy of the 
newly developed models. 


The Brouwer-Hori model is based on power series expansions 
in the satellite's radial distance from the center of the 
earth. It was found that these expansions do not converge. 
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Additional refinements in the drag force model should 
also be considered. The present theory considers the coeffi- 
oient of drag and ballistic number to be constants, but they 
may also be allowed to vary. The frontal area of a space craft 
in inertial hold or in a sun pointed orientation, could he 
modeled as a fourier expansion in the true longitude.- 

Finally, it is concluded that the analytical satellite 
theory need pot be limited by a simplified drag model. With 
the approaches developed in this report, it i? feasible to 
make the accuracy of the analytical theory competitive with 
preeisi-on numerical methods, while retaining a concise formu- 
lation and quickness of execution. 
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If one defines the fourier -series 


n 
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n (a-c aos i^ + as, sin I*) 

0 i=l 1 i 


m 


® (be 

j= 1 J 


cos j(p + bSj sin j(p) 


n+m 


C = Cg + (cCj^ cos k(f) + bs^ sin k(j)) 


where 


C = A*B and n>m 


then the coefficients of C are given by 


m 


c_ = 


=‘o’’’o + * 


Sc 

i=l 


ac . •be . + as . ‘bs . ) 

3. 1 1 


cc 


k 


^ •’’=k 


+ b 

0 


ac, 

k 




m 


s 

-; = 1 



0 



be. -as. ‘bs. 
J 1 3 


i+JT^k 

i+j=k 


+ 




0 

|i-J 

ac . ‘be ,+as . *bs 
L 1 J 1 jj 

|i“j i=k ' 





- 60 - 


cs, = 




■+ b„'0s. 


n 


r 

i=l 



0 

»l)s +as *bc, 
3 1 


i+j7^k 

i+J=k 


0 

— ^ ( ac . »bc . -as . ‘be . ) 

k '• 1 j X 2 ' 


|i-j l^k 
|i-Jhk 


Although the mathematical form of these equations appears 
to be quite clumsy, the algorithm can be programmed in a very 
concise and efficient manner. Therefore, the author has cho- 
sen also- to display the FORTRAN equivalent . 
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APPENDIX II 


As defined in reference 1 




(II. 1) 


or rewriting 


= 




_ e 


Q e oos (j) 


(II. 2) 


Therefore 


C? = ( q) (j) 


(II .3) 


The powers of a cosine can be written in- a fourier series with 

coefficients b!^ found by recurrence relations, 
k 


n 




oos <j) = / J B QOS k<f> 

k=0 ^ 


(II. 4) 


or 


n 


n n , 
e QOS (}) 


k=0 


„n ^n-k k , , 
e e QOS k(fi 


(II. 5) 


Since <j> is given by 


(}) = - (g+h) 
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0qTia^'iqn«v^ c^^.\be iwritrten . as 


n 


n 


oos 


n 


= L 


k-O 


B*" e“-'" 

n 


(C: 


OOS + 


s. 


.) 


sin ka ) (II. 6) 


where 

008 k(g+h) = .Sj 

6 k ~ ein k(g+h) = 




k-1 


(II. 7) 


and 


^ j = e oos <g+h) = Qp 2 
2 ^ (g+h) = -Qcr^ 


(II. 8) 


Therefore from equations (II. 6), (II. 3) and the recursive 
relations of (11,7) -one may write a complete expression for 
the powers of in jterms of a fourier expansion of the true 

longitude 

n 



